function fig2_18

global alpha lambda b

alpha=1;
lambda=1;
b=16;
t=115;
nx=500;
xmax=120;
x=linspace(-50,230,nx);
for ix=1:nx
    u0(ix)=exp(-x(ix)^2/b^2)*cos(lambda*x(ix));
    %h(ix)=int(@(s)exp(-s^2)*cos(s)*airy((x(ix)-s)/t^(1/3)),-5,5)/t^(1/3);
    h(ix)=(1/sqrt(2*pi))*rhom2(x(ix),t,lambda,b);
end

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.2,'MR',0.01,'ML',0.065,'P',0.04)

plot(x,u0,'b','LineWidth',1.1)
hold on
plot(x,h,'r','LineWidth',1.1)
axis([-40 xmax -1 1])

box on
grid on

set(gca,'ytick',[-1 -0.5 0 0.5 1])
set(gca,'yTickLabel',{'$-1$';'$-0.5$';'0';'0.5';'1'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

xlabel('$x$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('$h(x,t)$','Interpreter','latex','FontSize',14)
%legend('$\, t=0$','$\, t=40$','Location','NorthEast','FontSize',18,'Interpreter','latex','AutoUpdate','off')
legend('$\, t=0$','Location','NorthEast','FontSize',13,'Interpreter','latex','AutoUpdate','off')

%exportgraphics(gcf,'/Users/mark/Desktop/kg.pdf')

function RR=rhom2(x,t,lambda,b)
%  Integrate  over  [a, b]  using composite simpson and romberg
%  S(k) = composite simpson;  R(k) = romberg with composite simpson
%  err = relative iterative error using romberg

a=0; b=lambda+5/b;

%  interv = number of subinternals
tol=1e-8;
n=8;
S(1)=simp(a,b,n,x,t);
err=10*tol;
R(1)=1;  errR(1)=-1;
k=1;
while err>tol
    k=k+1;
    n=2*n;
    S(k)=simp(a,b,n,x,t);
    R(k)=(16*S(k)-S(k-1))/15;
    err=abs(1-R(k)/R(k-1));
end
RR=R(k);


function y=f(k,x,t)
global alpha lambda b
w=sqrt(1+k^2);
F =sqrt(2)*(exp(-(k - lambda)^2*b^2/4) + exp(-(k + lambda)^2*b^2/4))*b/4;
y=F*cos(w*t)*cos(k*x);


%%% composite simpson: n = number of subintervals
function s=simp(a,b,n,x,t)
xd=linspace(a,b,n+1);
h=xd(2)-xd(1);
s=f(xd(1),x,t);
for j=2:2:n
    ff2=f(xd(j+1),x,t);
    s=s+4*f(xd(j),x,t)+2*ff2;
end
s=h*(s-ff2)/3;














